R Notebook

Libraries

Load data

I am using the AHA dataset for writing my code because it is the simplest dataset I can work with. Simplest = significant number of baseline samples without issue of repeated measures, balanced distribution of classes, likely difference to be detected present.

a_ps = readRDS("/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/AHA/220714_MN00462_0226_A000H3WKM5/20221110_results/ps_objects/AHA_trnl_with_metadata_20221205.rds")

a_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 235 taxa and 95 samples ]
sample_data() Sample Data:       [ 95 samples by 33 sample variables ]
tax_table()   Taxonomy Table:    [ 235 taxa by 10 taxonomic ranks ]

c_ps = readRDS("/Users/aa370/Library/CloudStorage/Box-Box/project_davidlab/LAD_LAB_Personnel/Ammara_A/Projects/CHOICE/CHOICE_Trnl/CHOICE_20220912/ps_objects/choice_complete_20221205_clean.rds")

c_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 571 taxa and 171 samples ]
sample_data() Sample Data:       [ 171 samples by 30 sample variables ]
tax_table()   Taxonomy Table:    [ 571 taxa by 10 taxonomic ranks ]

Preprocess

Remove problematic samples

Taking out samples that are not listed to have completed the study and samples with more than 0 reads reduces samples from 95 down to 77.

Taking only the baseline timepoint brings further down to 44 samples

Keeping non-completed samples reduced the AUC of modelling substantially, remove these samples as before.

a_ps = a_ps %>%
  subset_samples(completion == "Completed" & reads >0) %>%
  subset_samples(timepoint == "Baseline")

a_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 235 taxa and 44 samples ]
sample_data() Sample Data:       [ 44 samples by 33 sample variables ]
tax_table()   Taxonomy Table:    [ 235 taxa by 10 taxonomic ranks ]

Class distribution in this dataset: Case = 20 Control = 24

sample_data(a_ps) %>%
  as.data.frame() %>%
  as_tibble() %>%
  group_by(group) %>%
  summarise(count = n())

Remove unassigned taxa

Debating whether this should be done for machine learning. Because we are simply throwing out because we don’t know what there are and even if they account for a small percentage of reads, they might still have important information.

I think I will refrain from removing unless I have a good reason to.

Don’t do this because some very important taxa for ML model turned out to be unassigned. - 12/8/22

Removing completely unassigned taxa reduces taxa from 235 to 82.

a_ps = a_ps %>%
  subset_taxa(is.na(superkingdom) == FALSE)

a_ps

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 82 taxa and 44 samples ]
sample_data() Sample Data:       [ 44 samples by 33 sample variables ]
tax_table()   Taxonomy Table:    [ 82 taxa by 10 taxonomic ranks ]

Merge phyloseq

combined = merge_phyloseq(c_ps, a_ps)

combined

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 674 taxa and 215 samples ]
sample_data() Sample Data:       [ 215 samples by 41 sample variables ]
tax_table()   Taxonomy Table:    [ 674 taxa by 10 taxonomic ranks ]

Aligning seq table

seqtab = combined@otu_table

seqtab_merged = dada2::collapseNoMismatch(seqtab)

combined = phyloseq(otu_table(seqtab_merged, taxa_are_rows=FALSE),
               sample_data(combined@sam_data), 
               tax_table(combined@tax_table))

combined

phyloseq-class experiment-level object
otu_table()   OTU Table:         [ 660 taxa and 215 samples ]
sample_data() Sample Data:       [ 215 samples by 41 sample variables ]
tax_table()   Taxonomy Table:    [ 660 taxa by 10 taxonomic ranks ]

sample_data(combined)

Sample Data:        [215 samples by 41 sample variables]:

Data transform

otu_clr = abundances(combined, "clr") %>% 
  t()

combined = phyloseq(otu_table(otu_clr, taxa_are_rows = FALSE),
                   sample_data(combined@sam_data),
                   tax_table(combined@tax_table))

rm(otu_clr)

Feature table

Extract data

Order of asv colnames in OTU table and tax table is the same.

features = combined@otu_table %>%
  as.data.frame() %>%
  as_tibble(rownames = NA) %>%
  rownames_to_column("sample")

feature_labels = c("sample")

features

Zero variance features

No features have zero variance

nzv = features %>%
  select(!sample) %>%
  nearZeroVar(saveMetrics = TRUE)

nzv %>%
  filter(nzv == TRUE)

rm(nzv)

Correlated predictors

summary(feature_cor[upper.tri(feature_cor)])

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.51377 -0.01487  0.03831  0.14088  0.18607  1.00000 

high_cor = findCorrelation(feature_cor, cutoff = .50) #using a more permissive cutoff because a lot of columns still left in the input table

filtered_features = features[,-high_cor]

rm(high_cor)
rm(feature_cor)

Linear dependencies

Input table

Modeling

Baseline performance

label_list = c("group", "group_shuffle")
results = loop_label(iterations = 20, input, label_list, percent = 0.8)

auc_df = results[[2]]

importance_df = results[[1]]

Loop

auc_list

 [1] 0.75 0.25 0.50 0.44 0.69 0.69 0.75 0.75 0.56 0.56 0.88 0.81 0.94
[14] 0.62 0.62 0.88 0.50 0.62 0.50 0.94

Plotting

AUC plot and t-test

Save external predictions

predictions

save = cummulative_external %>%
  #filter(auc == max) %>%
  left_join(samdf) %>%
  arrange(iteration, treatment, id, true_week)

Joining, by = "sample"

variable importances